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ABSTRACT 

"66": 

We introduce a simple and economical but effective method for including relativistic 
i— I ■ electron transport in multi-dimensional simulations of radio galaxies. The method 

Q_i' is designed to follow explicitly diffusive acceleration at shocks, and, in smooth flows, 

£^ . second-order Fermi acceleration, plus adiabatic and synchrotron losses for electrons 

in the energy range responsible for radio emission in these objects. We are able to 
follow both the spatial and energy (momentum) distributions of the electrons, so that 
direct synchrotron emission properties can be modeled in time dependent simulated 
7-H ' flows of this type for the first time. That feature is essential if simulations are to 

bridge successfully the fundamental physical gap between flow dynamics and observed 
Q\ ■ emissions. 

As an initial step towards that goal, we present results from some axis-symmetric 
Q\ . MHD simulations of Mach 20 light jet flows. These explicitly explore the effects of 

| shock acceleration, as well as adiabatic expansion and synchrotron aging in smooth 

O-f flows. The simulations demonstrate the importance of the fact that even for steady 

inflows jet terminal shocks are not simple, steady plane structures. Most importantly 
' this should play a very major role in determining the properties of synchrotron 

emission within the terminal hot spot and in the lobes generated by the jet back flow. 
In fact, the outflows are inherently complex, because of the basic "driven" character 
^ , of a jet flow. Consequently, the nonthermal electron population emerging from the jet 

may encounter a wide range of shock types and strengths, as well as magnetic field 
environments. 

We may expect to find a complex range in synchrotron spectral and brightness 
patterns associated with terminal hot spots and lobes. These include the possibility 
of steep spectral gradients (of either sign) within hot spots, the potential in lobes for 
islands of flat spectrum electrons within steeper spectral regions (or the reverse) and 
spectral gradients coming from the dynamical history of a given flow element rather 
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than from synchrotron aging of the embedded electrons. Finally, synchrotron "aging" 
in the lobes tends to proceed more slowly than one would estimate from regions of high 
emissivity. That is a consequence of the fact that those regions are ordinarily places 
where the magnetic fields are the strongest, so that the instantaneous rates of energy 
loss are atypical of the full history of the electron population. This feature supports 
earlier suggestions that nonuniform field structures may help to explain why dynamical 
ages of FRII sources often seem to be greater than the apparent age of the electrons 
radiating in the lobes, as measured in terms of spectral steepening, or absence thereof. 

Subject headings: particle acceleration - galaxies: jets - magnetohydrodynamics: MHD 
- radio continuum: galaxies 
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1. Introduction 

The standard paradigm for radio galaxies (RGs) is based on high speed plasma jets, formed 
in active galaxy nuclei (AGN), but penetrating far into circumgalactic environments and creating 
giant lobes of luminous material as a consequence of interactions with the intergalactic medium 



(e.g., Bridle 1992 ). The defining radio emission represents synchrotron radiation from relativistic 
electrons and magnetic fields. Those two crucial constituents are thought to be transported from 
the AGN, enhanced or generated by the jet-IGM encounter or both. Modern radio interferometry 



has provided richly detailed intensity, spectral and polarimetric images of RGs (e.g., Leahy e\ 



al. 1997| ), while rapidly advancing computational tools have allowed increasingly sophisticated 



multi-dimensional magnetohydro dynamical (MHD) simulations of the plasma flows (e.g., Clarkd 



1996). This paper presents the first such simulations that explicitly include time dependent 



transport of the relativistic electrons that produce the observed radio emission in such objects. 



Coleman & Bicknell 1988 carried out an early related calculation by mapping an analytic electron 



distribution onto a simulated bow-shock flow pattern, while |Jones fc Kang 1993| followed the time 



dependent behavior of a simplified electron distribution during the evolution of a shocked gas 



cloud. Clarke et al. 1989 and Matthews fc Scheuer 199C computed early models of synchrotron 



emission from RGs based on simulated fluid dynamical variables and intelligent guesses at the 
possible relativistic electron properties. But, no published simulations have followed the full 
electron distribution through a time-dependent evolution in a simulated RG-like flow. That 
important feature is a very difficult technical challenge, because the length and time scales needed 
to model effectively the microphysics of electron transport are typically many orders of magnitude 
smaller than those convenient to include in full scale hydrodynamical or magnetohydrodynamical 
simulations of RGs. 

To avoid these difficulties we introduce here a simple, economical method of relativistic 
electron transport that directly utilizes those mismatched scales and is suitable for time-dependent 
RG simulations. The method is designed to include the effects of diffusive shock acceleration and 
second-order Fermi acceleration, as well as adiabatic compression or expansion and synchrotron 
radiative cooling of the electrons in smooth flows. We present results of its initial application. 
There are many dynamical effects that are likely to influence the eventual synchrotron brightness 
distribution and spectrum to be expected in realistic models of RGs. Before it makes sense to 
attempt anything that could be termed a "real" model of a RG, it is crucial that the individual 
dynamical influences be understood. The present paper, therefore, is intended to begin addressing 
these individual factors. The role of shock acceleration of RG electrons is centrally important 
(e.g., [Heavens fe Meisenheimer 19871 ; piandford fc Eichler 1987| ; [Takarada 1989| ; jKriills 1992| ) 



While there are many studies of the physics of particle acceleration at plane or spherical 
shocks (e.g., |Kang &: Jones 199l|) , there are few previous simulations designed to explore shock 
acceleration in complex flows (but see e.g., Jones, Kang fc Tregillis 1994| ) and, as mentioned, none 



for flows similar to those expected in RGs. Therefore, our first priority is to explore how shock 
acceleration may be characterized and recognized in such environments, or the degree to which 
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such flows may tend to confuse simple interpretations of this process. 

Given these objectives and the fact that 2D, axis-symmetric flows are much easier to 
understand than 3D flows and still provide a first cut at the properties of RG flows, we limit 
this initial study to that problem. We have begun an extension of this work into fully 3D 
simulations and will present those results elsewhere. Since synchrotron cooling (commonly called 
"aging") is also likely to be a significant effect and will compete with and mask the influence of 
shock acceleration, we also include here an example of a flow that includes both of these effects 
together. We defer to later papers inclusion of second-order Fermi acceleration, not because it 
is necessarily unimportant, but because it is significantly more difficult to establish a unique, 
physical characterization of the momentum diffusion controlling it. Thus, at this stage, it would 
add significant complexity and uncertainty to an already complex set of questions. 

The plan of the paper is this. In §2 we outline our methods, including the new electron 
transport scheme. Section 3 introduces the parameters of the jet flows we have modeled, while 
§4 discusses our results. A brief summary of key findings is given in §5. An appendix contains a 
detailed discussion of practical requirements for numerical treatment of electron transport in RGs 
and justification for the scheme we introduce here. 



2. Methods 



2.1. Dynamics 



We evolve the equations of ideal nonrelativistic magnetohydro dynamics (MHD) in cylindrical 
coordinates (r,<f),z), with cj) an ignorable coordinate. All three components of velocity and 
magnetic fields are included, so the model is 2^-D. The code is an MHD extension of Harten's 



(Harten 1983) conservative, second-order finite difference "Total Variation Diminishing" scheme, 



as detailed in [Ryu fc Jones 1995j ; |Ryu, Jones fc Frank 1995| ; |Ryu, Yun fc Choe 1995| and [Ryu 
et al. 1998. The code preserves V • B = at each time step using an approach similar to the 



Constrained Transport (CT) scheme ( [Evans &: Hawley 1988| ) as described in detail by Ryu et 



al. 1998a. We use a passive "mass fraction" or "color" tracer, Cj, to distinguish material entering 
the grid through the jet orifice (Cj = 1), or from ambient plasma (Cj = 0). 



2.2. Electron Transport 

Our electron transport scheme is an adaptation of the standard diffusion-convection equation 
for charged particles (e.g., ^killing 1975 ); 

% = I»! (V •»-)-»■ V/ + V . («V/> + £| (/ D g) + Q, (2-1) 
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where f(x,p,t) is the isotropic part of the nonthermal electron distribution, k is the spatial 
diffusion coefficient, D is the momentum diffusion coefficient, Q is a source term that represents 
the net effects of "injection" and radiative losses at a given momentum, p, while the thermal 
plasma velocity is u. This equation is valid for "fast" (superthermal) particles when scattering 
is strong enough to keep the particle distribution almost isotropic and is appropriate on lengths 
large compared to the particle scattering length itself. 

The first term on the right accounts for adiabatic compression or rarefaction in the background 



flow. Eq. 1 2-1] is not valid within shocks, but when integrated across a velocity discontinuity, the 
first three terms account for "first order Fermi acceleration" at shocks, also known as "diffusive 
shock acceleration". The fourth term allows for second-order Fermi acceleration resulting from 



particle interaction with Alfvenic turbulence {e.g., |S killing 1975 ; Blandford &i Eichler 1987| ). 



As mentioned in the introduction and explained in detail in the Appendix, there is a severe 
mismatch between RG dynamical scales and diffusive transport scales that apply to eq. ][2-l|] for 
electron momenta relevant to radio synchrotron emission. That mismatch makes it impractical 



with conventional computational methods to solve eq. [2-1] as it stands to study electron transport 
in full RG flows. That is, the electron diffusive lengths and times are vastly smaller than those 
appropriate to the gasdynamics. The key physical fact is that the characteristic energy of electrons 
radiating synchrotron emission can be stated roughly as E ~ 3 v§ / B^ GeV, where vg is the 
observed frequency in GHz, and B\q is the magnetic field compared to 10 (j,G or nT. Thus, the 
electrons of interest mostly have energies <^10 GeV (p<^10 4 mc). The gyroradius of such particles 
is r g ~ 3 x 10 11 EceV / Bio cm - This leads to characteristic diffusion lengths and shock acceleration 
times in RGs <J AU and <^ yr, respectively compared to flow scales typically measured in kpc and 
kyr or larger. The most serious consequence of this comes from the requirement that numerical 
shocks must appear thinner than a particle diffusion length for accurate solutions of eq. |2-1| (see 
the Appendix). 

However, as also explained in the Appendix, that same mismatch can be exploited to develop 
a simplified equation of electron transport, provided f(p) is sufficiently broad that it can be 
represented as a piecewise power law over finite momentum bins. This last feature is very natural, 
in fact, in light of the same small time and length scales for electrons, since very rapid diffusive 
shock acceleration guarantees that GeV electrons emerge from shocks with power law momentum 
distributions on time scales much less than time steps required by the MHD. For much higher 
electron energies these simplifying conditions break down, but the method we derive could be used 
to provide a "low energy injection spectrum" for those particles, as well. 

For continuity with the remainder of our discussion we present here the simplified transport 
equation, but refer readers to the Appendix for derivation details and tests of its validity. We divide 
the momentum range of interest into a modest number, N, of logarithmically spaced bins bounded 
by p , . . . ,p N . Defining m = lnpi/p , bin widths can be expressed as, Ay, = y i+1 - yi = lnp i+ i/pi. 
We can integrate equation ][2-l[| within each momentum bin to define = 4tt JJ l+1 p 3 fdlnp as the 
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number of electrons in the bin. For convenience we normalize rij by the total plasma mass density, 
to form bi = rii/p, then assume a piecewise power law f(p) = fa (p/pi)~ qi within pi < p < Pi+i, 
so that rtj is given in terms of fa, % and Ay^ by eq. [ |A2| |. For these relatively low energies 
diffusive transport at shock discontinuities is properly handled by defining the form of the electron 
distribution just down stream of a shock to be the steady state power law appropriate to the jump 
conditions for that shock. That fixes the ratio bi+i/bi at shocks. Electrons are injected from the 
thermal plasma at shocks using a common injection model, so that a fixed fraction, e, of the total 
electron flux through a shock is injected and accelerated to the appropriate power law momentum 
distribution. Away from shocks, in smooth flows, but where the diffusion lengths of the electrons 
are much smaller than dynamical lengths, spatial diffusion is negligible. We can easily include the 



effects of synchrotron "aging". When all of these features are included eq. [2-l| becomes (see also 



eq. jSH) 



dbi (\ qD 1 p\p 3 f Pi+i 



r_i = 4vr[-Vn-^ + — -\— , (2-2) 
at \3 p z T so pJ p Pi 

where t so = | %£jj g g defines the synchrotron cooling time at momentum p, which we take 

arbitrarily to be p = 10 4 mc. In this expression, ax is the Thomson cross section, and 

Ub = B 2 /(8tt). Note that the right hand side of eq. \2-2\ is the difference between fluxes at the 



two momentum boundaries of the bin i. Thus, we have a conservative, finite volume scheme that 
depends on a simple model of sub-grid structure (in momentum space) for its accuracy. 

Although straightforward to include, we now drop the second order Fermi term, "D" from 
these initial simulations for the following reasons. First, in the vicinity of shocks, this acceleration 
process is likely to be much less efficient than first-order Fermi acceleration. That is simply 
because the second-order process comes from nearly balanced energy fluxes of waves propagating 
parallel and antiparallel to the magnetic field that resonantly scatter with the electrons; i.e., on 
"isotropic" Alfvenic turbulence. First-order Fermi acceleration, on the other hand, depends 
only on waves propagating in the streaming direction that will be generated by the streaming 
particles themselves. In either situation resonant waves have wavelengths comparable to the 
particle gyroradius. Downstream of the shocks where particle streaming may be less important, 
it is significant for such low particle energies that the resonant wavelengths are probably within a 
couple orders of magnitude of thermal ion gyroradii (taking u s ~ 10 8 cm s _1 ). Such waves may 
be fairly strongly dissipated, for example, by nonlinear Landau damping (e.g., [Lee fc Volk 1973 ). 



Thus the level of relevant isotropic Alfvenic turbulence will likely depend on a local cascade of 
hydrodynamical turbulence (e.g., Grappin et al. 19821 ). While we will certainly want to understand 



what role post-shock, second-order acceleration might play (e.g., Borovsky fc Eilek 1986 ; Kriill 



1992[ ), only ad hoc models could be applied at present. Since this set of computations represents 
the first effort to treat nonthermal electron transport inside time dependent models of RGs, and 
since there are many possible competing influences to be understood before "realistic" models are 
sensible, the best initial strategy is to restrict ourselves to the most straightforward physics that is 
likely to be important. First-order acceleration is efficient and almost certain to happen at shocks, 
so its role must be understood from the beginning. 
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2.3. Computation of Synchrotron Emissivity 

From the spatial distribution of bi and qi along with the MHD variables, p and B, it is 
straightforward to compute the synchrotron spectral emissivity, j u . Ignoring corrections due to 
slow curvature in the electron spectrum, j u , can be expressed as (cf. |Jones, O'Dell fc Stein 1974) 

ju = jao^ff( P )p q {^ff^ (2-3) 

with a = j ao ~ 1 (cf. Pones, O'Dell fc Stein 1974| ) and vb ± = vbcosO the electron cyclotron 
frequency in terms of the magnetic field projected over the angle 6 onto the plane of the sky. 
f(p) along with q are found from mapping bi and g, onto the critical synchrotron frequency as 
p = mcJ • We- have interpolated q between momentum bin centers to give a a continuous 



form. Eq. [2-3] would break down for strongly curved spectra, but by that point in calculations 
of the type used here the local index will be too steep for the emission to be significant. In a 
calculation designed to predict the surface brightness distribution of a RG model one should 
include explicitly the effects of variations in cos 9. That is not our intention here, and such a 
calculation would need to evaluate carefully the effects of "sub-grid" fluctuations in the field 
directions and of the orientation of the jet with respect to the line of sight. Our aim here is 
to understand how electron acceleration and transport are likely to be reflected in the overall 
emission properties of a region within a RG, independent of the location of the observer. For 
those purposes it is better to ignore the factor cos#, so will in this paper simply take it to be 
unity. Again, we are not in this paper computing surface brightness distributions from these 
simulations, which would involve line-of-sight integrations of j v . That would be inappropriate, 
given the axis-symmetric nature of the simulations, and would further exaggerate the dynamical 
constraints imposed by this symmetry. 



3. The Simulated Jet Properties 

Our initial simulated MHD jets are all dynamically identical. They have a simple, "top hat" 
velocity profile, are Mach 20 with respect to the uniform ambient medium (Mj = Uj/c a = 20), 
in thermal pressure balance with it, and have a density contrast rj = Pj/p a = 0.1. The jet enters 
at z = 0, with an initial radius of 36 zones, while the entire uniform grid is 384 x 1536 zones 
(10|rj x 42|rj). Defining length and time in units of initial jet radius (rj = 1) and ambient sound 
speed (c a = y/jP a / 1 Pa = 1, with 7 = |), the simulations are followed for about 11.7 time units, 
when the bow shock of the jet reaches the right z boundary. Reflecting boundaries are used along 
the jet axis, while continuous boundaries are used elsewhere. There is a background poloidal 
magnetic field {Bbackground) (Bfj, = B r = 0; B z = B zo ), with a magnetic pressure 1% the gas 
pressure; i.e., = iP, with j3 = 10 2 . The in-flowing jet carries an additional toroidal magnetic 
field component appropriate to a uniform axial jet current with a return current on the surface of 
the jet; i.e., B^ = 2 x B zo (r/rj) for r < rj. At r = rj,/3 = 20. The in-flowing jets are slightly 
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over-pressured at the outside. 

These simulations are truly MHD rather than HD, despite the apparent weakness of the 
initial magnetic field. Even at the initial strengths modeled here the magnetic fields play a 
role in the evolution, especially in the cocoon. There are numerous locations where the plasma 
(3 ~ 10, so that its pressure is not entirely ignorable. There is, however, a more important and 
much less recognized role for the field in flows initialized with even (3 » 1. In 2-D and 3-D 
MHD simulations of the Kelvin-Helmholtz instability, for example, we have demonstrated crucial 



dynamical contributions from magnetic fields in flows with initial (3 > 10 3 (e.g., Jones et al. 1997 



Ryu et al. 1998b| ; |Jones, Ryu fc Frank 1998| ). The point is that through field "stretching", 



magnetic tension can be increased in complex flows to the point where it is comparable to or 
exceeds the Reynolds stresses of the gasdynamical flow ( coming from spatial variations in puiUj), 
even though the magnetic pressure may be smaller than the gas pressure. The relative importance 
of Maxwell to Reynolds stresses can be roughly represented by the Alfven Mach number in a 
flow, Ma = u/va- When the local Ma drops to small values in a complex flow, that can lead 
through reconnection to so-called "dynamical alignment" and flow "self-organization" . There are 
also many locations within our model flows where M^^l. In short, the flows are made smoother 
by the presence of the field, even though in the original configuration the field nominally appeared 
too weak to influence the dynamics. A comment may be in order on the meaning of magnetic 
reconnection obtained from an "ideal MHD" code. While reconnection is most fundamentally 
a topological transition (e.g., , Axford 1984), it cannot take place without resistive, dissipative 
effects on small scales. In an ideal MHD code these dissipative effects are numerical in origin, so we 
cannot model the microphysics of driven reconnection (which is still not understood at this time, 
in any case). However, there is evidence that the phenomenology of reconnection is often correctly 
captured. For example, Miniati et al. (1998), using the Cartesian version of the code we employ 
here, demonstrate clear development in 2-D supersonic MHD cloud simulations of the so-called 
"resistive tearing-mode" instability that is physically associated with driven reconnection. 

We present here three examples of the electron transport within the dynamics of the above 
jet flows. Their properties are listed in Table |l[ Electrons are modeled in the momentum range 
Po = m e c and pn ~ 1.63 x 10 5 m e c. All three include the effects of diffusive shock acceleration. 
Models 1 and 2 are alike in that the electrons have negligible influence from synchrotron aging 
, while Models 1 and 3 are alike in that the electron populations originate entirely with the 
in-flowing jet. Model 2 differs in that the electron population is mostly injected locally from 
thermal plasma at shocks within the modeled flow. Model 3 differs, then, in that electrons are 
influenced very significantly by synchrotron aging. 

In the "adiabatic" Models 1 and 2, there is very little curvature to the electron spectra, so 
in those models we used four momentum bins (i.e., N = 4) giving In = 3. (See the Appendix 
for justification of these numbers.) In Model 2 we injected electrons at po using the simple model 
described in §A.3. Electrons should be injected at slightly suprathermal postshock energies, so the 
value po = m e c was selected as appropriate for jets with speed Uj ~ 0.1c, assuming that shocked 



Table 1. Summary of Simulations 



Model a N b In-flowing Electrons Jet Spectrum d Shock Injection Cooling Time f Bbackground! 

(bi) q je t (a j e t) (e) (t so ) ^Gauss 

1 4 10~ 4 4.4 (0.7) 0.0 40 10 

2 4 10~ 6 4.4 (0.7) 10~ 4 40 10 

3 8 10~ 4 4.4 (0.7) 0.0 4 30 



a All models used Mach 20 jets (M, = Uj/c a — 20), with a density contrast, ij = pj/p a = 0.1. Units come from 
a jet radius, rj = 1, an ambient density, p a = 1, and a background sound speed, c a — •yPa/pa = 1 (7 = 5/3). 
The cylindrical computational box was 10.67 Tj units high and 42.67 Tj long. End time for each simulation was 
t = 10.67 units. There was a uniform, poloidal jet magnetic field, B z o that also extended to the background, so 
B z o = B background with P — P/Pb = 100, while the incoming jet also contained a toroidal field, — 2 x B z o(r/rj). 

b Number of electron logarithmic momentum bins spanning the range po = 1 mc to pjv = 1.63 x 10 5 mc. 
(lnpjv/p = 12) 

c Ratio of nonthermal to thermal electrons in the incident jet flow. 

d qjet is the momentum distribution power law index of the in-flowing electron population; aj et is the corresponding 
synchrotron index. 

c Assumed fraction of the thermal electron population injected to the nonthermal population as it passes through 
a shock. 

f Time for electrons to cool below momentum, p = 10 4 mc in the ambient, or background, magnetic field, measured 
in time units rj/c a . 

g Field that would produce the listed cooling time, assuming illustrative values for the simulated jet: Uj = 0.1c, 
Vj — 1 kpc. This would also map the displayed synchrotron emissivity data onto 1.4 GHz. Applying the additional 
constraints on the jet listed above, the implied jet density and kinetic luminosity are pj ^ 3 x 10~ 29 Bl ackground g 
cm~ 3 and K kinetic = f p^rf w 1.2 x 10 43 Bl ackground erg sec -1 , with B background in ^Gauss. 
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electrons are thermalized to a temperature comparable to but smaller than the ions. Since the 
nonthermal, cosmic-ray electrons are passive, all results could be scaled for different choices of p$. 

To include synchrotron aging it is necessary to relate the characteristic cooling time t so 
defined in association with eq. [A3| to the computational time units, here given by rj/c a . The 



cooling time depends as well on the magnetic field strength and the momentum p as 

t s = 2.5 x 10 1 — — (3-1) 



Va Mj r jk B 



2 ' 
10 



where t s = r so 2 p 4 i s the electron momentum in units 10 4 mc, ii,-g is the jet speed expressed in 
units of 10 s cm s _1 , Mj is the jet Mach number (here set to Mj = 20), rjk is the incoming jet 
radius in kpc, and Bio 1S the strength of the magnetic field in units of 10 /iG, or equivalently in 
units nT. By comparison, the end points of our simulations are t = 10|. 

For Models 1 and 2 = 40, the objective being to produce negligible aging 

for electrons of interest during the simulation. To project emission from electrons with momenta 
p^W 4 mc into the radio band, we should be considering magnetic fields J>10 fiG. Assuming, for 
example, that B z q = 10 /iG, and Ujg = 30, then rjk ~ 1-1, so that the length of our computational 
box would be about 45 kpc, and the time unit, rj/c a » 7 x 10 5 years. 

By contrast, for Model 3 we set r so = 4.0, so that synchrotron aging is now very significant 
for electrons in the energy range of interest over the time scale of the simulation. Keeping the 
same physical jet speed and radius, as well as the time unit, this case would correspond to an 
axial magnetic field B z $ ~ 30 /j,G. Other associated dependent parameters are outlined in Table 
1. These combinations are illustrative only and not meant to be particularly representative of any 
real RGs. Other combinations would work as well for our present purposes, which are aimed at 
the physics of the electron transport. For the simulation of Model 3, we used eight momentum 
bins (N = 8) to cover the same range as the other models, as justified in §A.4. 



4. Discussion 

4.1. Flow Dynamics 

We first outline some key dynamical properties of the simulations, remembering that the 
MHD properties of the different electron transport models are identical except for dimensional 
scaling factors. Fig. 1 shows at t = 10.67 images of the plasma density, p, the magnetic field 
intensity, B, and gas compression, = — V • u. The compression highlights shock structures. 
These images remind us that such flows are complex and unsteady. All of the structures are 
ephemeral or at least highly variable. Jets, as strongly driven flows are not reasonable to model 



as equilibrium structures (cf. Sato et al. 1996). As described by many before us {e.g., Norman, 



[Smarr &: Winkler 1985| ; |Lind et al. 1989 ), once equilibrium is broken, jet plasma alternately 



expands and then refocuses while interacting with its self-generated cocoon, creating oblique jet 
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shocks as a result. ^] Those shocks are neither steady nor stationary, however. Most important to 
our discussion, the oblique jet shocks interact episodically with the terminal jet shock, causing 
it to vary in both strength and structure. The jet never re-establishes an equilibrium in its flow. 
The terminal shock also includes a "Mach stem", so that in 2-D at least, jet flow coming down 
the outside of the jet, near rj, always exits through an oblique shock. At times there is little or no 
perpendicular terminal shock; then, most of the emerging jet flow is only weakly shocked. When a 
perpendicular portion to the terminal shock exists, it usually is strong, with a compression ratio 
a ~ 3.8 (— » q ~ 4.08; a ~ 0.54) in these simulations. The many other shocks and the oblique 
segments of the terminal shock are mostly weak (Fig. lc), but they can exhibit density jumps, 
a £ 2 - 3 (-> q « 6 - §). Along the jet axis intersecting oblique shocks can be strong, but these 
intersections intercept a small fraction of the flow, so will not produce a large population of flat 
spectrum electrons. 

Vorticity is shed from the outer edge of the terminal shock to form the turbulent jet cocoon. 
There are distinct episodes of strong "vortex shedding" coincident with disruption and reformation 
of the terminal shock. Once shed, the vortices interact with the Kelvin-Helmholtz unstable 
boundary layer of the jet, generating an even more complex back-flow and perturbing the jet flow, 
as well. All of these dynamical features are represented in the Fig. 1. snapshot. The large "rolls" 
visible in Fig. 1 are remnants of vortex rings that were shed earlier than the time displayed. The 
magnetic field in the back-flow is dominated within an axis-symmetric flow by B^, both because 
this magnetic component is enhanced by stretching in the expanding flow emerging from the jet, 
and because magnetic reconnection of the B r , B z field components annihilates most of that field 
from the back-flow in an axis-symmetric geometry. This flux annihilation inside vortical flows is 
well-known and sometimes termed "flux expulsion" in the MHD literature ( [Weiss 1966 ). In a 



3-D flow both the vortex and magnetic field structures would be stretched, twisted and tangled; 
.e., the back-flow would become turbulent and disordered on small scales (e.g., Norman 1996| ; 



i 



Ryu et al. 1998b; Jones, Ryu & Frank 1998) leading to vortex and magnetic flux tube complexes. 



4.2. Electron Transport &: Emissivity 

The electron transport we model takes place "on top of" this dynamics; i.e., , the electrons 
do not feed back dynamically on the flow. To check this for consistency we confirmed that for the 
numbers we assume, the electrons never represent more than about 0.1 % of the kinetic energy 
density in the jet or the back flow. Remember that the dynamics behind the electron transport 
and synchrotron simulations discussed below is identical in each model. The electron distribution, 
and especially its momentum dependence, reflect the flow history of a given fluid element as 
appropriate for the particular transport assumptions, as well as the current fluid condition. Our 



5 In this paper "oblique" and "perpendicular" refer to the angle between the shock face and the velocity field, not 
the shock normal and the magnetic field, as is customary in the particle acceleration literature. 
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objective in this paper is to begin to learn how to interpret the emergent patterns. Since the 
astronomical tool for this is synchrotron emission, we will explore the physics in that paradigm 
using local emissivities. We expect that this approach will make much more tractable the 
interpretation of full 3D simulations. For the convenience of the discussion we express emissivities 
as they would be at 1.4 GHz, based on the fiducial magnetic field values listed at the end of §3 
and in Table |]. They can be easily rescaled by choosing other combinations of source parameters. 
Illustrative examples of the emissivity properties in these models are shown in Fig. 2-4. Fig. 5 
explores the correlations between emissivity, spectral index and magnetic field strength using the 
same data displayed in the other figures. At the start we point out the obvious "limb brightening" 
of the emissivities seen in all of the models displayed. This feature is not characteristic of real radio 
lobe surface brightness in FRII sources, but is an artifact of the assumed axis-symmetry of the 
simulations; particularly the persistence of the large "rolls" in the back-flow and the decoupling 
that occurs in axis-symmetry between Bj, and (B r , B z ). In an equivalent 3-D simulation the 
bright structures should be better mixed into the turbulent interior regions of the back-flow. This 
is one reason we defer any display of model surface brightness until we have 3-D structures in the 
models. That does not detract from the present discussion, however, since this initial exploration 
is mostly aimed at identifying physical relationships in emission properties and flow properties, as 
well as recognizing emissivity characteristics that we may eventually be able to use as diagnostics 
of the history of the local electron population. 



4.2.1. "Adiabatic" Models 1 and 2 

Fig. 2a and 2b provide grayscale images of \ogj v at 1.4 GHz for the "adiabatic" electron 
models 1 and 2 at t = 10.67. For these two models emissivity values from eq. R2-3H are displayed 
when they are within a factor 3 x 10 3 of the peak emissivity. Smaller emissivities are blacked out. 
For comparison we also display with the same dynamic range in Fig. 2c a "pseudo emissivity", 

3 

based on the MHD properties of the flow; i.e., j c = CjPB^, modeled after the approach 
introduced in Clarke et al. 1989. We have added the jet color variable, Cj to Clarke's original 



definition, because of differences in our assumed background magnetic field. Clarke's magnetic 
field vanished outside the jet flow and its cocoon, whereas we used a simpler model that continues 
the axial field into the background. Since Cj = in the background our definition allows j c > 
only for regions containing jet plasma, so that they are comparable to Clarke's. On the other 
hand, for the Model 1, 2 and 3 emissivities based explicitly on transported relativistic electrons, 
the absence of a significant electron population in the background makes this factor immaterial. 
Fig. 3a and 3b show with an inverse gray scale the 1.4 GHz spectral index distributions, 
a = — 31ogjy/(?log i/|i.4 = (q — 3)/2, {a\ and 0:2) ° ver the same emissivity range in Models 1 
and 2. 

Fig. 2 compares emissivities produced by the two "adiabatic" models represented by j\ and 
j2- to the "pseudo emissivity", j c . The three emissivity images have clear similarities. All show 
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a dominant "hot spot" , down stream (to the right) of the jet terminal shock. Each produces 
qualitatively similar patterns of bright emission within the back-flow. The jet is illuminated in 
the j c model emissivity with patterns that resemble those of the j\ model. Of course, the jet is 
largely invisible in j'2, because the jet electron population is assumed very small in that model. 
The model congruencies reflect the strong roles of two key ingredients to the emissivity; namely, 
the strength of the magnetic field and adiabatic compression and expansion. Note, however, 
that ji and ji have much greater local contrast than j c . Since the magnetic field distributions 
are identical for all three emissivity calculations, differences in the effective radiating electron 
distributions must be responsible. That is apparent by noting within smooth flows that m oc p » , 
so ji 2 oc niB^~ oc p~ . This is actually rather similar to j c oc p ' B ' , except for variations 
modeled in bi = rii/p and qi, but not modeled by j c . The formula given for j c is based on an 
assumed a = 0.5, so contrast could be increased in j c by assuming a significantly steeper index. 
The natural extension of j c to general a would be j c oc p 2a p l ~ 2a (B sinip) a+1 v~ a (D. Clarke, 
private communication). For a = 1 this gives j c oc (PB) 2 /p, just as for ji 2 with a = 1. That 
would not represent the same effect noticeable in j\ and j'2, however, since from Figs. 2, 3 and 5, 
we can see that large ji 2 come primarily from regions where a < 0.7. 



4.2.2. "Synchrotron Aged" Model 3 

In Fig. 4 we show the analogous simulated synchrotron emissivity and spectral information 
for our Model 3, which includes the effects of synchrotron aging . Recall that, except for this 
influence, Model 3 is identical to Model 1. Thus, Fig. 4a (j'3) can be compared to Fig. 2a 
(ji), while Fig. 4b (03) can be compared to Fig. 3a (cci). In addition, we illustrate the level of 
spectral curvature of Model 3 in Fig. 4c, using of the difference between the simulated spectral 
indices at 1.4 GHz and 5.0 GHz (63 = a^{f>GHz) — a^lAGHz)). Note that we use formal 
definitions of a rather than the traditional observed value defined by the ratio of fluxes measured 
at two distinct frequencies. The intrinsic "observed" index can be found by the simple formula 
a\ A = a + 0.5 * 5s- These "observed" indices would be almost the same as a near the origin of the 
jet and in flows just either side of the terminal shock in Fig. 4. On the other hand the "observed" 
spectra would tend to be ~ 0.05 steeper in mid regions of the jet and in most of the back flow, 
cocoon regions bright enough to show in this figure. But, we mention another caveat with regard 
to comparisons with directly observed properties. Beyond other complications already mention, 
including differences between emissivity and surface brightness, a real observation would likely 
smooth over a physically significant range of sight lines. Non-black regions in Fig. 4 cover a wider 
range of emissivity values (namely a factor 3 x 10 4 ) than in Figs. 3 and 4, since the inclusion 
of synchrotron aging adds considerable contrast to the emissivity distribution. This factor 10 
extension captures approximately the same physical regions as for the other models. 

The emissivity distribution, j'3, for Model 3 including synchrotron aging, is shown in Fig. 
4a. As mentioned earlier, the dynamical range found in j'3 is greater than for any of the other 
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models, since there is considerable steepening visible in the spectral index as one moves away from 
electron sources; namely, the jet origin and strong shocks. Steepening due to aging is apparent 
even within the jet, where 03 increases from 0.7 at the origin to values near 1 even before the first 
oblique shocks are encountered. Although the emissivity distribution within the jet in Fig. 4a is 
very similar in appearance to that for Model 1 in Fig. 2a, the factor 10 increase in dynamic range 
covered in Fig. 2a slightly de-emphasizes the stronger emissivity decreases caused by synchrotron 
aging. By contrast to the adiabatic models, where the synchrotron spectral index remains near 
a ~ 0.7 almost everywhere in the jet, the oblique jet shocks do produce a noticeable flattening 
in the spectra of Model 3. The oblique shocks have compression ratios, cr<^3 corresponding to 
a postshock index, q > 4.5; a > 0.75, so they only flatten an incident spectrum that has been 
steepened by radiative losses. In fact, from Fig. 4c it is apparent that in Model 3 there is 
significant spectral curvature in the synchrotron spectrum even within the jet. The nature of that 
is also visible in the electron momentum distribution plot in Fig. 7c. While the oblique jet shocks 
are relatively ineffective as electron accelerators, places where these oblique shocks intersect form 
strong shocks with their normals aligned to the jet axis. They produce spectra synchrotron spectra 
as flat as a ~ 0.52 (see Fig. 5). But, as already mentioned, the portion of the jet flow intercepted 
by these shocks is very small, so in these simulations they have little impact on the large scale 
synchrotron emission within these flows. 



4-2.3. General Comments - Electron Transport & RG Dynamics 

We begin this section with some observations about patterns in the synchrotron emissivity 
spectra within the jet structures themselves. In Fig. 3 we see that the jet spectra for the adiabatic 
Model 1 and 2 are almost everywhere a = 0.7, which was the injected spectrum. That is despite 
the existence of several jet shocks and apparent "knots" in the flow before the terminal shock, 
where the spectrum flattens to a ss 0.54. Should we have expected to see spectral flattening at the 
jet "knots" ? In these two models the answer is no, since shocks only flatten an incident electron 
spectrum if it is steeper than the index q s = 3a /(a — 1), where a is the compression ratio through 
the shock (see §§2.2 and A. 3. 2). Thus, only shocks with compression ratios, a > 3.14 would lead to 
flattening within these jets. In fact, except very near the jet axis the jet shocks are oblique, with 
compression ratios, o < 3, in these simulations. There is one very restricted flatter spectrum region 
in the jet in Models 1 and 2 right on the axis and just to the left of the terminal shock (Fig. 3). 
There a very strong shock has formed at the intersection of oblique shocks. However, as mentioned 
earlier, the volume of jet plasma involved is very small, so this feature is hardly noticeable. On the 
other hand, synchrotron cooling in Model 3 steepens the incident electron spectra coming into 
the "knots", so that in Fig. 4b we can see clear evidence for spectral flattening at the jet "knots", 
albeit to indices, aj>0.7, in these simulations. The details of such features will probably vary with 
the details of the assumed jet structures and Mach numbers. For example, if the jet flows and 
magnetic fields conspired to produce the greatest emissivities close to the jet axis, then spectral 
changes at shocks would be more apparent, since jet shocks do tend to be transverse on axis. 
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The observational evidence b caring on this is still limited. Our choice of q Q — 

4.4 (a = 0.7) 

was intended to represent a "typical" FR II RG jet spectrum, of course. It is most often hard 



to isolate structures in RGs that are clearly "only" jet flows (see, e.g., Rudnick & Katz-Stone 



1996 ) because jets are often relatively weak emitters and because they can, therefore, be confused 
by contamination from the cocoon surrounding them. This makes it especially hard to isolate 
changes from the jets to the knots, because the contamination from cocoon emission (expected to 
be relatively steep-spectrum) is greater for the inter-knot regions than in the knots. In some cases 
where the jet data look convincing, such as for M87 (actually a FR I source), the radio spectrum 



seems to be virtually constant along the entire jet with no indication of curvature (e.g., Biretta 



et al. 1991; Meisenheimer, Roser fc Schlotelburg 199"6| ). The extension of this spectrum to the 



visual band is what strongly constrains aging in M87, since the estimated cooling times for radio 
emitting electrons are longer than typical estimates of propagation times. And, indeed, there 
are indications, despite the constancy of a in the radio band, that the radio-to-optical spectrum 
steepens between the "knots" in the M87 jet ( [Sparks, Biretta Macchetto 1996[ ). In Cyg A the 
jet radio spectrum appears also to be locally nearly straight, but to exhibit different slopes at 
different locations ( Rudnick fc Katz-Stone 1996 ; Katz-Stone et al. 1993j ). Other radio jets may 



show evidence of steepening of the spectrum along the jet (e.g., [Mack et al. 1998|) , but data do 
not yet allow detailed analysis in terms of electron transport. Clarke et al. 1986 reported evidence 
for spectral flattening at a knot in the jet of Cen A. The spectral index flattens to a ~ 0.7, from 
the surrounding emission where a ~ 1.5. If these changes really are within the jet flow and not 
influenced by contamination, this case could easily be accounted for by the effects we see in our 
simulations. 

The spectral index distributions in Fig. 3 show for our adiabatic models that much of the 
bright emission outside the jet itself, especially in the hot spot, is associated with spectra flatter 
than the a = 0.7 expected from the in- flowing jet plasma. Even the brightest parts of the cocoon 
are dominated by spectral indices with a^0.7. That must be largely due to acceleration in the 
terminal shock, of course. Closer examination reveals considerable fine structure, however, in the 
spectra. Fig. 5 shows, especially for Model 2, that reasonably bright emission is produced with 
spectral indices ranging roughly over a wide range, a ~ 0.54 — 0.75. The a ~ 0.54 emission comes 
from electrons accelerated in the strong, plane portions of the terminal shock. Near the head of the 
jet in this model the brightest emission in the back-flow still tends to have relatively flat spectra, 
a <^ 0.6. However, it is interesting to note that the spectra of the brighter areas in Models 1 and 
2 are often steeper farther back in the cocoon. That does not come from radiative losses, which 
are negligible in these two models, but from increased mixing between flat and steep electron 
populations. This is an important effect to recognize, since it mimics what is usually assumed 
to be due to synchrotron aging in the back flow. The specific locations of flat and steep spectra 
in these examples, should be viewed with caution, of course, since once again we are observing 
axis-symmetric flows that do not mix in the same ways as the more general 3-D flows will do. 



Where do these steep-spectrum electrons come from in the adiabatic models? The steeper 
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spectrum emission in the hot spot and the cocoon represents plasma that has escaped passage 
through the strong, plane terminal shock. One such path is plainly visible for Model 1 in Fig. 
3a. There we can trace emission with a = 0.7 all the way from the jet origin on the left into 
the cocoon back- flow. There are even a = 0.7 "fingers" of unmodified jet electron populations 
just outside the jet and near the left of the grid. The source of these electrons is made apparent 
by the invisibly low j% values (Fig. 2b) in the same regions. That plasma flows into the cocoon 
through the outer, oblique portion of the terminal shock where the shock compression a < 3, so 
the spectrum is unmodified from its original form in the absence of synchrotron aging. We note 
in an axis-symmetric flow that plasma emerging from the jet terminus near the outer radius of 
the jet must initially flow back into the cocoon closer to the jet channel than plasma emerging 
along the jet axis, since 2D flow "streamlines" do not cross. That also explains why near the hot 
spot flatter spectrum emission tends to extend along the outer, forward parts of the cocoon in 
Models 1 and 2, since that represents plasma passing through a strong terminal shock nearer the 
jet axis. As noted, turbulence within the back-flow eventually mixes these populations together, 
steepening the average spectrum. There is also a related very important source of steep spectrum 
electrons in the cocoon in these adiabatic models; i.e., there are time periods with little or no 
perpendicular terminal shock. Then most or all of the emergent jet plasma carries electron 
populations unflattened by the terminal shock. This is especially relevant when there is fresh 
electron injection at shocks, as in Model 2. There, the steeper spectrum emissivities can be 
introduced by electron injection in oblique portions of the terminal shock. 

Finally, in the context of the adiabatic models we mention that there is significant, bright 
emission in Model 2 coming from electrons in the shocked, ambient IGM; i.e., where Cj = 0, 
so j c = 0. This comes from IGM electrons injected to the relativistic population by the jet bow 
shock and now embedded in regions with moderately strong magnetic fields. Fig. 3b shows that 
the spectral range is large in these regions, with a ^ 0.6. 

At first glance the emissivity distribution for the synchrotron aged Model 3 in Fig. 4a is 
very similar to the adiabatic models, especially Model 1, which shares the same electron injection 
history. Even accounting for the greater dynamic range displayed in this is especially true 
near the terminal hot spot. Within the hot spot that similarity extends to the spectra, which 
are relatively flat in both models. But, closer examination reveals some interesting differences 
that come from the effects of strong synchrotron aging. First, even within the hot spot, Model 3 
displays a substantial spectral range, from close to a ~ 0.5 to values greater than a ~ 1. In fact, 
the brightest regions of the hot spot have spectral indices just greater than a = 0.7, by contrast to 
values near 0.54 for Models 1 and 2 (see Fig. 5). The strong spectral steepening that takes place 
at the very front of the terminal hot spot region is a product of the very strong magnetic field that 
forms at the nose of the emerging jet. This is evident in Fig. 4e and 4f, where it is apparent that 
the emission from the strongest fields (which lie in this region), while very bright, is also quite 
steep. From Fig. 5 we can see that the strongest fields in the hot spot are ~ 8 times greater than 



Bbackground, so from eq. |3-lj , the cooling times of the emitting electrons are quite short; namely, 
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r s ~6x I0~ 2 rj/c a . In that situation, with the absence of any additional acceleration within the 
postshock flow, spectral steepening is a virtual certainty. In fact, strong spectral index gradients 
in some RG hot spots have been identified using "spectral tomography" techniques (Treichel & 



Rudnick 1998j), showing that real hot spots have structures that may eventually enable us to 



explore electron radiative histories and origins, once we understand the behaviors of more realistic 
models. 

For the most part the brighter regions in the cocoon of Model 3 have spectra with aj>l, 
and on the whole, the cocoon has a steeper spectrum than the hot spot. These properties are in 
rough, general agreement with observational expectations (e.g., Dennett-Thorpe et al. 1997j ; Mack 



et al. 1998 ). It is interesting to note that there are "pockets" of flat spectrum emission in the 
back flow, however. There is a particularly prominent region in the forward-most swirl just above 
the terminal hot spot, visible in Fig 5. Electrons in these regions were subjected to strong shock 
acceleration, but were shed in a vortex event without ever passing through a region of strong 
magnetic field capable of causing strong aging. The field in their current position is actually below 
^background ( see Fig- lb), so their cooling time has become a large fraction of the full simulation 
time. These regions are moderately bright because the electron density is relatively high. As 



possible evidence that such structures might form in real RGs we mention that [Mack et al. 199S 



point to the existence of "flat spectrum islands" (q ~ 0.6) in the lobes of NGC 315 where a ~ 1.3. 

The presence of flatter spectrum populations within the back flow is also related to the fact 
that additional steepening of the spectrum of the cocoon is slight in the brightest regions in 
Model 3 once one looks past the most forward portions (see Fig. 4b). From Fig. lb and Fig. 5f 
it is plain that the magnetic fields in the visible regions farther back in the cocoon are roughly 
comparable to those in the forward cocoon or the much of the terminal hot spot. In particular 
they are greater than Bf }ac j tgroun ,i. Thus, the nominal cooling times based on current conditions 
would be t s < Ar,j/c a , while the age of the jet at this time is t ~ 10.67rj /c a . The reason for this 
apparent contradiction between short cooling times and slow aging is the highly unsteady and 
nonuniform nature of the cocoon flow, especially the magnetic field, in the cocoon. Because of that 
the electrons spend only a fraction of their lives in strong field regions where they emit strongly. 
Thus, they cool much more slowly than one would naively predict, based on emissivities. This 
supports previous suggestions that the radiative ages of RGs appear less than their dynamical 
ages because of small magnetic field filling factors (e.g., Myers fc Spangler 1985| ) or filamentary 



magnetic field structures (e.g., Eilek et al. 1997 ). 



We can recognize some additional and generally interesting behaviors from Fig. 5. For 
example, from Fig. 5a and 5c for most of the points representing shocked jet emission (including 
all the points with a < 0.7 in Fig. 5a) there is an envelop restricting the strongest emission to 
the flattest spectra. Note, for example, in Model 1 none of the brightest emission outside the jet 
comes from regions with a > 0.6, despite the fact that the magnetic fields in the fainter, steep 
spectrum emitting regions are as strong as in the bright, flat spectrum regions (see Fig. 5b). 
The "nose" of the (j u ,a) distribution in Fig. 5a, 5c, 5e corresponds to the terminal hot spot 
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emission, where fields are strong and electron densities are high. The largest emissivities naturally 
correspond to the strongest field regions. The converse is not true, however; strong fields do not 
necessarily produce a large emissivity. The other ingredient in j v , of course, is the electron density, 
f(p). Fig. 5 is, therefore, showing us a wide range of electron densities associated with magnetic 
field strengths contributing to the emission, and even some degree of anti-correlation between field 
and rii in the back-flow. Back-flowing plasma has mostly undergone adiabatic expansion after 
emerging from the jet, while the component of the magnetic field has been enhanced in some 
locations by stretching of field lines at the same time. In Model 3 (represented in Fig. 5e,f), 
where synchrotron aging also acts to reduce emission by steepening the spectrum, the patterns are 
actually very similar to those in the adiabatic models. However, much of the emission is shifted 
upward to steeper spectra; roughly by Aa ~ 0.2. The points remaining near a ~ 0.5 correspond 
to shocked points within the jet near its axis, and points just past the jet terminal shock, where 
aging effects are not yet important. Finally, we comment that more generally, there are many 
"micro patterns" visible in Fig. 5. Those tend to trace hydro dynamical features, so that they are 
revealing histories of the electrons. Unfortunately, those flow patterns are usually unsteady so 
that it would be difficult, if not impossible to use them to determine explicit connections between 
individual points and flow dynamics. 

5. Conclusions 

We have developed a simple but effective numerical scheme for transport of low energy 
relativistic electrons suitable for studies in time dependent simulations of radio galaxies. The 
scheme can follow the evolution of electrons accelerated by first order Fermi acceleration in shocks, 
second order Fermi acceleration, adiabatic cooling and synchrotron cooling, a.k.a. "aging", in 
smooth flows. From the electron spatial and momentum distribution and the magnetic field 
properties found through the fluid evolution one can directly compute the radio synchrotron 
emission. In this paper we have applied those methods to axis-symmetric MHD jet flows including 
first order Fermi shock acceleration and synchrotron aging effects on the electrons, in order to 
identify some of the behaviors that may be expected in multi-dimensional flows. We begin with 
2D flows since they are much simpler than 3D flows to study and to understand. We have focussed 
our discussion on effects that should occur in both 2D and 3D flows, so that what we learn here 
will facilitate future studies with more general flows. Briefly, the most important findings of this 
work that we believe to be independent of our specific simulation assumptions are: 

• The fact that jet terminal shocks are not simple, steady plane structures should play a 
major role in determining the properties of synchrotron emission within the terminal hot spot 
and in the lobes generated by the jet back flow. In fact, the outflows are inherently complex, 
because of the basic "driven" character of a jet flow. As pointed out recently by [Sato et al. 1996 , 
Complexity is basic to all driven non-equilibrium plasma flows. RGs are dramatic examples of 
such flows, so they should be influenced by episodic patterns of instability and reorganization. In 
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particular, the strength and form of the terminal shock is constantly and sometimes dramatically 
changing. Thus, the nonthermal electron population emerging from the jet may encounter a wide 
range of shock types and strengths, as well as magnetic field environments. 

• Because of the previously mentioned behaviors, we may expect to find a complex range in 
synchrotron spectral and brightness patterns associated with terminal hot spots and lobes. These 
include the possibility of steep spectral gradients (of either sign) within hot spots, the potential 
in lobes for islands of flat spectrum electrons within steeper spectral regions (or the reverse) and 
spectral gradients that result from the dynamical history of a given flow element rather than from 
synchrotron aging of the embedded electrons. 

• Synchrotron "aging" in the lobes tends to proceed more slowly than one would estimate 
from regions of high emissivity. That is a consequence of the fact that those regions are generally 
places where the magnetic fields are the strongest, so that the instantaneous rates of energy loss 
are atypical of the entire history of the electron population. This feature is apparent even in 
axis-symmetric flows; it should be stronger in real, 3D flows, since for equivalent initial structures, 
local magnetic field enhancements should be greater in 3D (Jones, Ryu fc Frank 1998|) . O ur 



finding supports earlier suggestions that nonuniform field structures may help to explain why 
dynamical ages of FRII sources often seem to be greater than the apparent age of the electrons 



radiating in the lobes (e.g., Eilek et al. 1997 ) 



There are many uncertainties inherent in attempts to model such complex objects as RGs. 
The flows are inherently intricate and many key parameters are largely unknown at present. That 
is why we have examined first the simplest systems that may be qualitatively representative and 
have tried to restrict our questions to those not wholly dependent on uncertainties coming from 
the limitations imposed by our various choices. Even from these few initial results we can see the 
importance of modeling electron acceleration and transport and the significance of shock history 
to the radio spectra in RG hot spots and lobes. The interpretation of synchrotron spectra in 
RGs clearly is not simple, but represents the combined influence of several dynamically coupled 
factors. That should, perhaps point us to another analysis approach; namely, thinking directly 
in terms of "Complexity" itself. (Sato et al. 1996| argue that in such systems we need to shift 
our paradigm from a focus on specific elements (e.g., the existence of a terminal shock that can 
accelerate particles) to the consequences of interactions among elements (e.g., long range order 
and the relationships between dynamical structures). The task then becomes one of understanding 
how the "microphysics" (e.g., particle transport) can help to establish relationships among the 
various system elements. 
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A. Appendix 

A.l. Constraints on Computational Diffusive Electron Transport in RGs 

A number of studies have considered time dependent solutions to eq [|2-1|] for charged particles 



in shocked flows (e.g., 


Falle & Giddings 1987; Drury 1991; Kang & Jones 1991; Kang, Jones & Ryu 


1992 




Berezhko, Ksenofontov & Yelshin 1995; 


Kang & Jones 1995; 


Kang & Jones 1997 


). Those 



studies and most others like them have focussed on the transport of nonthermal ions. Although 
diffusive transport is distinct for nonrelativistic electrons (e.g., Levinson 1996| ), relativistic 



electrons can be treated exactly like ions of similar rigidity, so the above treatments are directly 
relevant to the transport of electrons responsible for synchrotron emission. All of the papers 
mentioned, however, considered 1-D flows for a single shock system. This limited application 
reflects a severe numerical constraint in solving eq. [2-1], especially for low energy particles. With 



existing numerical schemes it is not practical to carry out a full and direct solution to eq. [2-l| 
within a "large-scale", complex multi-dimensional flow. The reasons are straightforward. First, 
note that eq. j2-l|/ applies only outside of shocks, since its validity is restricted to particles whose 
scattering lengths and gyroradii are large compared to the shock thickness, which should be a few 



thermal ion gyroradii (e.g., Kang Jones 1995 ; Kang &: Jones 1997 and references therein). Its 
application to diffusive shock acceleration depends entirely on being able to match upstream and 
downstream solutions properly at the shock (e.g., Drury 1983; Drury 1991). An accurate numerical 



solution to eq. [2-1] requires that the shock be a discontinuity or close to it. Numerical schemes 
for compressible fluid dynamics, including the one we employ, spread shock structures over at least 
one or two zones. One of the necessary matching conditions for diffusive electron transport is that 



the distribution / be continuous across the discontinuous physical shock. So, if we apply eq. [2-1 



at a shock, we must demand, to avoid errors within the (unphysical) numerical shock, that the 



advective flux term in eq. [2-1] integrate to the same value as at a true discontinuity. It is simple 
to show that the errors in this are first order in Ax s /xd, where Ax s > Ax is the numerical shock 
thickness, Ax is a zone width at the shock and xj = k/u s is the diffusion length for energetic 
particles at the shock, with u s being the velocity jump across the shock. Therefore, we must have 
Ax s /x d « 1. 

To see what constraint this gives in a RG application we need to estimate Xd- There is 
empirical evidence that near shocks the diffusion coefficient, k, is a steeply increasing function of 
momentum with a form resembling the so-called Bohm diffusion coefficient, k ~ r g w, where w 
is the particle speed (e.g., Ellison et al. 1993). For GeV electrons in a magnetic field exceeding 
~ 10 /iG, for example, we would estimate for Bohm diffusion near a strong shock with speed, 
10 8 cm s , that k<^10 cm s and x^lO cm. Lower energy electrons, especially those 



-21 - 



near shock injection energy, would give smaller Xd by at least two orders of magnitude. Thus, a 
direct simulation in a RG of the full eq. [ 2—1 1 for electrons even at GeV energies would require 
numerical shocks much thinner than 10 14 cm. To include electrons down to injection energies 
in this situation we would necessarily require numerical shocks to be no thicker than about an 
AU. This is a spatial constraint, so temporal sub-cycling techniques to deal with associated rapid 
acceleration are not sufficient to deal with this problem. Since the simulation scale in a RG is 
typically tens or hundreds of kpc one would need to resolve shocks on a scale no more than about 
10~ 6 of the main grid scale and much smaller than this to capture the particle behavior with 
eq. [ 2-1 1 near injection energies. For simple, 1-D flows this kind of refinement near shocks may 
be feasible using a nonuniform grid, although one still must deal with numerical issues coming 
from the accompanying orders of magnitude range in the time steps required for the solution. An 
adaptive mesh refinement (AMR) approach may eventually make the analogous multi-dimensional 
treatment feasible, but until those codes are available for MHD, including diffusive transport, an 
alternative way around the problem is necessary for meaningful calculations. 



A. 2. A Practical Approach to Dealing with Electron Transport 



Fortunately, for the problem at hand we can actually exploit the small diffusion length 
constraint to introduce a simplified version of eq. [2-1] that still correctly captures the essential 
physics of electron transport, at least when the dynamical feedback of the electrons can be ignored; 
that is when the electrons can be seen as "test particles." [Jones fc Kang 1993 applied a very simple 



version of such a scheme several years ago, and |Jun fc Jones 1998| have briefly described a variant 
of the one discussed below. To begin, we remind readers of the point in §A.l that we should 
distinguish between behaviors at shocks and within smooth flows, away from shocks. Therefore, 
we divide discussion of our proposed scheme into these topics and begin with shocks, since the 
emergent properties of the electron distribution help to determine our approach in smooth flows. 



A.2.1. Transport Through Shocks 

At shocks we must properly balance the change in advection at the shock against diffusion 
away from it. From the arguments in the preceding section, we must find a way to make the 
numerical shock much thinner than a conventional numerical zone width. The most obvious 
numerical solution to this problem is to treat the shock as a real discontinuity; i.e., confine 
the shock transition to a zone boundary as far as particle acceleration is concerned. Generally, 
as mentioned earlier, that involves imposing the matching conditions to / across the shock 



{e.g., Drury 1985; Drury 1991) to properly account for both diffusion and advection on both 



sides of the shock. Drury 199l| used this method to find analytic time dependent test particle 



solutions for / at a plane shock, for example. Berezhko, Ksenofontov fc Yelshin 1995| used it 



numerically to treat time dependent ion acceleration at the spherical blast wave of a supernova 
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remnant, including backreaction on the local flow. Those particular methods are not practical to 
apply to numerical schemes in general flow patterns, because they are computationally intensive 
and require highly customized assumptions. However, for the low energy electrons in RGs that 
are of interest to us in this paper, it is possible to accomplish the required outcome in a manner 
that is easily adapted to general numerical schemes. The key points are: 1) as already discussed, 
diffusion lengths are very small compared to typical numerical zones on either side of the shock 
for the particles under consideration; i.e., Ax >> Xd] and 2) the characteristic diffusion timescale, 
td = x d/ u s adjacent to the shock is much smaller than numerical time steps associated with the 
standard CFL condition; i.e., At ~ Ax/u s >> td- In fact, At/td ~ Ax/xd, so, from our earlier 
discussion, both conditions are easily satisfied by many orders of magnitude in RG simulations. 

The first point means that the entire upstream "precursor" to the shock that controls 
evolution of / is well-contained within a single numerical zone. This prohibits treatment of 
backreaction by the nonthermal particles on the flow, since that depends on knowing the velocity 
gradient within the precursor. But, electrons accelerated in shocks are generally thought not to 
produce significant backreaction (e.g., Blandford & Eichler 1987). In the work we present here 



the nonthermal electron population is not dynamically coupled back to the flow. The second 
point also means that the precursor is formed instantaneously by comparison to the numerical, 
dynamical time step. More directly, it means that the timescale for particles to be accelerated 
at the shock to the energies of interest is also very small compared to the numerical time step. 
Assuming for simplicity that k is the same upstream and downstream of a strong shock, the 
mean time for particles to be accelerated from an initial momentum, p Q , to momentum, p, can be 
written approximately as t a « ^ JJ o n(p')d\n.p' ( |Lagage fc Cersarsky 1983[) . For Bohm diffusion 
of relativistic electrons near a strong shock t a ~ 20 x td ~ 2 x W 7 E GeV /(B 10 u 2 s8 ) sec, where 
u s s is the shock speed expressed in units of 10 8 cm s _1 . For GeV electrons this is about a year 
for our selected characteristic numbers, compared to a typical CFL-determined dynamical time 
step in a RG simulation J>10 3 years. The synchrotron radiative cooling, or "aging" timescale, 

t s ~ 10 15 /(f 9 2 Bi ) sec, should be comfortably longer than either of these for the particles of 
immediate focus, so we can deal with effects of radiative cooling separately on large scales. It has 
been demonstrated analytically ( Drury 1991| ) and numerically (Kang & Jones 1991) that after a 



time At a particle distribution f Q flowing into a steady plane shock emerges downstream with the 
appropriate steady state solution for p/p a t « 1, where p a t is the momentum that can be reached 
by diffusive acceleration during the time interval At, as defined by t a . For our application, where 
p a t(At) is much larger than momenta of interest if At is a typical CFL hydrodynamical time step, 
it becomes unnecessary to treat the detailed time evolution of / at shocks with eq. p-l| . Instead 
for GeV electrons we can assume immediately downstream of a shock that in response to diffusive 
shock acceleration / oc p~ q , where q = min(q s ,q ), with q Q = — din f Q /dlnp the logarithmic slope 
of the in-flowing electron distribution at momentum p, and q s = 3r/(r — 1), the standard diffusive 
shock acceleration power law slope. In this r = P2/P1, the compression ratio through the shock 
(e.g., Prury 1983Q . 
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The above discussion also leads to the conclusion that, at least as it emerges from shocks, the 
electron momentum distribution for energies of direct relevance for radio emission will generally 
be well-approximated by a power law. That emergent form may change in subsequent propagation 
over finite lengths due to radiative losses, mixing between particle populations with different initial 
spectra or second-order Fermi acceleration. Nonetheless, this suggests an obvious simplification to 
the energy distribution that greatly reduces the effort required to solve eq. [2-1] outside of shocks 
in RG flows. In particular q = — din f /dlnp should be a slowly varying function of \np. Then, 
for numerical purposes it makes sense to approximate the momentum distribution as a piece-wise 
power law over bins of finite width in Inp. The practical coarseness allowed for such a momentum 
grid will be determined by the expected curvature of the spectrum and the information needed 
from the spectrum in the momentum range of particular interest. 



A. 2. 2. Transport in Smooth RG Flows 

In smooth flows, eq. [ [2-1(1 is directly applicable. Here the relative importance of advective and 
diffusive fluxes is given by a ratio u5x/k, where u is a characteristic advection speed across zone 
boundaries and 5x is the scale length of V/ in the smooth flow. Assuming 5x > Ax, it is easy to 
show with the diffusive properties mentioned earlier, that on a typical numerical grid with ~ 10 3 
zones spanning the RG dimensions in any direction, u5x/k > uAx/k >> 1. Thus, advection of / 
for low energy particles in the smooth flows should always dominate over spatial diffusion even 
if the spatial diffusion coefficient is several orders of magnitude greater than that we discussed 
next to the shocks. We conclude that we can neglect the spatial diffusion term in eq. |2-1|/ for low 
energy electrons in RGs except at shocks. 



A. 3. A Simple Electron Transport Scheme for RGs 



To solve eq. [2-1 we should divide a momentum range of interest into N bins bounded 
by po, . . . ,pn- It is most practical to work with logarithmic intervals in p, so we generally use 
yi = In pi/ po, giving bin widths, Ayi = yj +1 — y% = \npi+\/pi. For a conventional, finite difference 
treatment of eq. [2-1] it is necessary to use bins with Ayi « 1 (e.g., Kang fc Jones 1991 ). 
However, with the conservative, finite volume approach defined here, it is practical to use much 
greater bin widths, Ay^l, when the smooth behavior of \nf(p) outlined in §A.2 applies. 



A. 3.1. Smooth Flows 



To begin, we can integrate equation [2-1 within each momentum bin to define 
rjj = 47r fp l+1 p s fdlnp as the number of electrons in the bin. It is convenient (but not 



necessary) to normalize rii by the total plasma mass density, to form b{ = rii/p. Then eq. [2-1 



-24- 



becomes, for one spatial dimension the conservative, finite volume equation 

dbi „ fldu D a In A p 3 f P<+ i d (. .dbi\ Q' 

47T - — + ' ■ 



\3dx p 2 dy J p Pi dx \ K ^dx) p ' 
where < k >= ^ f ^f d P qi _ ^ f Pi + 1 p 2 Qdp, and d Idt is the Lagrangian time derivative. No 

J p^\7fdp J Pi 



approximations have been added to obtain eq. [|A1| from eq. 2-1 1. 



Now we explicitly include the two simplifying features into eq. |A1| that enable us to treat 
low energy electron transport efficiently in RGs. We shall also define Q' specifically, so that it 
accounts for synchrotron losses and shock injection of low energy electrons. Remember that eq. 
p-l| or fAlfl applies only to the smooth flows. In those regions we assume the piecewise power law 
f(p) = fi (p/Pi)~ qi within pi <p< p i+ i, so that 

n l = ^M-[l - (-«-)«-3]. (A2) 
Qi ~ 3 Pi+x 

Since f(p) is continuous we can also write fi + \ = fi{pi/pi+i) qi ■ That enables us to find the fluxes 



at momentum boundaries in eq. | A.1 to account for adiabatic expansion, second-order Fermi 



acceleration or radiative cooling. Then also recall from §A.2 that we can neglect the spatial 



diffusive transport term in eq. |A1] for low energy electrons within smooth portions of flows for 



the conditions expected in RG simulations. 

Synchrotron aging of electrons is straightforward to include in this formulation when there is 



rapid angular redistribution of electrons (a necessary condition for the validity of eq. [2-1]). Then 
at momentum p there is a flux A7rp 2 p s f, with p s = —^—p, where t so = f ™\j B \ defines the cooling 
time at some convenient momentum p (see, e.g., pones, O'Dcll Sz Stein 1974 ). In this expression, 



ot is the Thomson cross section, and Ub = B /(8tt). To account for inverse Compton cooling 
from the cosmic microwave background B 2 should be modified to B 2 + B 2 , where B^ = 3.2 pG at 
the current epoch. Putting these features together, we are left in smooth flows to solve the simple 
transport equation 

djH =A fldu_ q D + ^p_\p*l 
dt \3 ox p l T so p) p Pi 



Since eq. [A3] is in conservation form, our main constraint is the accuracy of the fluxes 
written on the right hand side. Given an initial set of 6, and q^, along with underlying plasma 
properties, those can be found accurately at the beginning of the time step, within the piecewise 
power law approximation by using eq. [ |A2[[ and the continuity of /. The fluxes can easily be time 
centered, by using the slope of /, to achieve second order time accuracy in the update. One must, 
of course, include a constraint on the time step, At, beyond the MHD CFL condition. That can 
be expressed conservatively to account for adiabatic, synchrotron and second-order acceleration 
effects by the combined condition, At < max (Ax/a, T so p / p^ ,p 2 / (qD) min ) , where a is the fastest 
MHD signal speed. The first condition, for adiabatic expansion, corresponds to the standard CFL 
condition. The other two need to be examined explicitly, but so long as the synchrotron cooling 
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and second-order acceleration take place on "hydro dynamical timescales" , they will also usually 
be enforced by the CFL condition. 

In general qi will also change during a time step. That set of values is updated from the new 
distribution 6j by inverting eq. [ |A2| , and again using the continuity of /. An appropriate set of 
boundary conditions is also required, of course. We expect spectral curvature will be common 
at the highest momenta, but usually minimal at the lowest momenta. So, for the simulations 
presented here we have assumed that qi is continuous at po, but that the spectral curvature is 
continuous at p^. The solution to qi must be done iteratively. As pointed out by pun fc Jones] 



1998| , it is possible in principle to find qi exactly if the momentum bins are of uniform size and one 
applies the boundary condition q Q = q±. For significantly curved spectra, however, we have found a 
simple, approximate method to be relatively fast and sufficiently accurate in these applications. In 
this method we assume for the inversion of eq. |A2[ to find that q is a linearly varying function 
of y between the centers of adjacent momentum bins; that is, there is uniform spectral curvature 
between adjacent bin centers. Thus, qi is interpreted as the average q (or midpoint value) within 
the interval Then it is straightforward to design a Newton-Raphson iteration scheme 

to solve for each qi. In practice we use only a single iteration from the initial guess 

ln%ti 

qi (trial) = 3-—^, (A4) 

which would be exact for a pure power law on uniform bins. 

The number of bins necessary in using eq. [A3| over a given momentum range will depend 



on the amount of spectral curvature expected. Jun & Jones 1998 showed if one can assume that 



f(p) is a pure power law at any given physical location, then two momentum bins are sufficient to 
follow the evolution of / accurately, accounting for adiabatic expansion, advection and consequent 
mixing of different power laws. In our application that is not adequate, since we may expect 
significant deviations to develop from pure power law form. On the other hand, from thermal 
injection energies to those ~ GeV the relevant A hip < 20, while a range \Aq\ = 2Aa = 3 would 
cover electron populations producing synchrotron spectra with slopes 0.5 < a < 2, assuming 
q > 4. Greater curvature would involve spectral slopes, q > 7, so that very few particles would 
be represented at the higher energies. As demonstrated below we have found that it is possible 
to capture the evolution of electron spectra utilizing only a few momentum bins for practical 
problems. The computational cost of following each electron momentum bin is roughly comparable 
to following each dynamical variable in the flow (of which there are eight in MHD). The same 



statement would apply for a conventional scheme to solve eq. [2-1], such as that we used in Kang 



fc Jones 1991c but there, a much larger number of momentum bins would be necessary and much 
finer spatial resolution near shocks would be needed to obtain meaningful solutions in RGs (see 
§A.l and comments at the end of §A.4). 
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A.3.2. Shocks 

To model "diffusive shock acceleration", following the arguments of §A.l, we impose at 
detected shocks the analytic solution outlined in §A.2. That is, the emerging distribution has a 
form f(p) oc p~ q , where q = 3r/(r — 1), with r the compression ratio through the shock provided 
q s < q, where q is the in-flowing value in a momentum bin. Normalization of / is determined 
by the total number of electrons; that being the sum of the incident population in the interval 
[PO; Pn] and any injected electrons as defined below. We ignore shocks with density jumps r < | 
(q s > 7; a > 2), since they contribute very few nonthermal particles. 



For fresh injection at shocks we adopt a commonly used simple model (e.g., Kang & Jones 



1991); namely, we inject a small, fixed fraction, e, of the thermal electron flux through the shock. 
Thus inside a shock we solve the equation db/dt = Q' in j/p, where b = J2^i and 

Qinj = • (A5) 

W s is the Lagrangian shock velocity, [i e is the electron mean molecular weight, mn is the mass of 
the proton. As for electrons flowing in from upstream, shock injected electrons are distributed 
in momentum as a power law with the appropriate slope from diffusive acceleration theory, as 
described before. 



A.4. Tests 

We have carried out numerous tests of the method described above and find it to be both 
very stable and to produce consistent results, provided roughly lAg^+il < 0.75. In the simulations 
described here we have constrained jA^^il < 0.67. Our various tests included 1-D adiabatic 
advection of analytically defined convex and concave spectra, as well as shock tube simulations. 
We carried out a variety of 2-D simulations, such as spherical blast waves and obliquely intersecting 
shocks, where we could compare analytic expectations of the resulting electron spectra. We also 
tested the ability of the scheme to handle synchrotron aging by following the evolution of an initial 



power law electron distribution compared to the well-known analytic result of |Kardashev 1962 
(but, most commonly associated with |Jaffe fe Perola 1973 ) for synchrotron aging with rapid pitch 
angle redistribution. 

A 2-D application of this scheme to the propagation of a young SNR impacting an interstellar 
cloud is presented in Jun & Jones (1998), and can serve to demonstrate its ability to capture 
expected behaviors correctly for a spherical blast wave, as well as for mixed flows. We offer 
here another simple 2-D test calculation that could be used by other workers for quantitative 
comparison. It involves the interaction between two obliquely intersecting plane hydro dynamical 
shocks, and is illustrated in Fig. 6. Moderate strength shocks are more challenging to capture 
correctly with regard to electron acceleration, so the example is designed to create several shocks 
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with Mach numbers in the range 1.5 — 10. There are five dynamical regions, as outlined below 
and as shown in Fig. 6a. Region 1 is at rest, with p\ = 1 and P\ = 0.6, giving unit sound speed. 
Region 2 (also uniform) lies behind a Mach 10 right-facing shock initially placed vertically in the 
middle of the box. The initial conditions in region 2 are, thus, p2 = 3.884, P2 = 74.88, u X 2 = 7.425, 
u y 2 = 0. Region 5 similarly lies behind a Mach 4 shock, which, in this case is propagating 
diagonally up and to the left. This shock was initially placed so that it bisected both boundaries, 
and, thus, was just touching the vertical shock. Initial conditions in region 5 were, /O5 = 3.368, 
P 5 = 11.85, u x5 = -1.989 and u y5 = 1.989. 

Propagation of each of the two original shocks into the postshock flows of the other creates 
new regions 3 and 4, which are separated by a slip line. We carried out this test on a 512 x 512 
Cartesian grid. Except for the right edge of the box standard continuous boundaries were used. 
On the right these were modified to include the symmetry of the flow there, thus preserving its 
self-similarity. The whole flow pattern is self-similar, in fact, except near the bottom of the box, 
where it is not possible to construct suitable boundary conditions. That leads to the curvature in 
the shocks visible there, and the incipient vortex at the bottom of the slip line. Using shock polars 
(e.g., , Courant & Friedrichs (1976)) it is simple to compute that regions 2 and 3 are separated by 
a Mach 1.51 shock, while regions 4 and 5 are divided by a shock with Mach number 3.66. 

For this particular test there was an initial cosmic-ray electron population that was a uniform 
1% of the thermal electron population in each region. The cosmic-rays arbitrarily were given an 
initial, uniform power-law form with index, q = 4.5, corresponding to original injection at a Mach 
3 shock. Many other spectral choices could be made, as well. There was no injection of cosmic-ray 
electrons at the shocks during the simulated flow from this point, although it is straightforward to 
include that effect, as well. 

Fig. 6b illustrates the computed distribution of spectral indices after the flow has evolved. 
Since the flow is self-similar, except for the mentioned boundary effects, the chosen time does 
not matter, so long as it is sufficient to allow significant material to pass through each shock. 
As mentioned the electrons initially had the same spectrum everywhere. Particle acceleration at 
the shocks during their propagation during the simulation has modified the electron spectra in 
regions 2', 3", 4' and 5'. The steady-state spectrum behind a Mach 1.51 shock (region 3) from a 
monoenergetic incident flux would be q = 7.11. Since that is steeper than the initial spectrum 
there is no change between regions 2 and 3'. However, region 3" was first processed by the Mach 
10 shock, so we find, as expected, electrons there to have the same form as those in region 2'. 
From the analytical jump conditions we predict the spectral indices in regions 2', 4', and 5' to be 
respectively, qy = 4.040, q& = 4.323, q^ = 4.267. In this numerical example they come out to 
be q2> = 4.041, q^ = 4.333, and q^ = 4.299, with variations less than or of order 10~ 4 . Thus, 
the errors in q, which are naturally greatest for oblique shocks on a Cartesian grid, are all less 
than 1%. Since there is little mixing and no synchrotron aging during this test calculation, all 
the spectra should very good power laws, which we confirmed using 4 momentum bins covering 
Amp = 20. 



-28- 



This test also demonstrates the ability of our routine to advect the electrons cleanly in 
the different momentum bins, since the spectral indices are computed from the actual particle 
distributions among the distinct momentum bins. So, for example, the population differences 
between regions 3' and 3" in the higher momentum bins are more than an order of magnitude, 
and that difference is maintained accurately with a transition over about three spatial zones. Any 
error or significant diffusion there would lead to a spurious flattening of the spectrum along the 
boundary region 3'. The spectra there computed for the two highest momentum bins is almost 
precisely the same as for the two lowest bins. Similar comments apply to the slip line between 
regions 3" and 4'. Recall, of course, that we assume for these calculations that real, physical 
diffusion is very small on scales of the computational zones. 

To demonstrate further the behavior of our electron transport scheme we include Fig. 7, 
which shows a snapshot of the results from three 2D test simulations exactly like the Model 3 jet 
described in §3, except that the spatial resolution here was only half as fine. In these calculations 
an electron population with po =mc, pn ~ 10 5 mc and a power law index q = 4.4 enters with the 
jet on the left of the grid. That population is modified by adiabatic and synchrotron cooling and 
by any shocks encountered along the way. The grayscale image shows the log of the total electron 
number near the end of the simulation. For numerical convenience, the computation also includes 
a small population of steep spectrum {q = 7) background electrons in the ambient medium. They 
are faintly visible in the image, especially through their compression at the jet bow shock. Fig. 7 
also shows the computed electron momentum spectral index at four selected locations on the grid. 
Each plot gives the values of qi found in simulations with the number of momentum bins, N set 
to 4 (triangles), 8 (stars) and 12 (squares). 

Note first the plot (b) in Fig. 7, which shows the spectrum of the background population, 
changed only in response to synchrotron aging in a uniform background magnetic field. The solid 
curve gives the prediction, q(p,t) = —dlnf/dlnp from Kardashev 1962| ; namely, 



t v 

g(p,*) = go + (go-4) (A6) 

TsoP 



where qo is the initial slope. In this case qo = 7, and the time shown was t = t so . Equation 
A6 is valid, of course, only for < 1, since higher momentum particles have all cooled below 



TsoP 



Pcool = P^-f 2 -- In this case p coo i = p = 10 4 mc. All three numerical results agree with eq. [A6] in the 



range of momenta where the curvature is less than the numerical imposed constraint mentioned 
earlier; namely that the change in % between adjacent bins is less than |. For N = 4 violation 
occurs above p ~ 10 3 mc in this case, but for both N = 8 and N = 12 the numerical results are 
in reasonable agreement below about ^p C ooi- The deviation at higher momenta has little practical 
consequence, since the numerical slope there is still so steep and the associated 6, so small that 
it will have little influence on the evolution at lower momenta (vanishing flux contribution) and 
would contribute negligibly to emission processes. 



Similarly with each of the other three particle spectra shown (all corresponding to the 
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dynamical time displayed in the image), the agreement between the N = 8 and N = 12 solutions 
is very good until the computed spectra are too steep to influence observable properties. For 
the two spectra (c) and (d, where curvature is small below p ~ 10 4 mc, the N = 4 solution 
agrees quite well with the other two. We conclude that the scheme is effective and reasonably 
robust. When there is only modest spectral curvature in the momentum range of interest, as 
few as 4 momentum bins may be adequate to capture the evolution of the spectrum. For the 
first two RG models presented in §3, synchrotron aging is negligible, so we use iV = 4 in those 
calculations. In situations where strong curvature needs to be captured, such as those involving 
strong synchrotron aging , a larger number of bins are needed. In the Model 3 presented in this 
paper we consider emission almost entirely from electrons having p < 10 4 mc, and we judge there 
that N = 8 is adequate for the issues under discussion. This ability to capture basic behaviors 
of electron transport in RG flows with only a few momentum zones results from the conservative 
nature of the scheme, the quasi-power law properties of the electron distribution and the fact that 
the microphysics of the electron transport can be analytically well-approximated in smooth flows 
and at shocks. The impact of this fortuitous outcome is clear when one considers that any scheme 
to follow the electrons will require a computational effort per variable tracked that is comparable 
to the effort necessary to follow a dynamical variable, such as the gas density. In this scheme we 
need 2N variables (pi, with N ~ a few, so that the net effort to follow the relativistic electrons 
is comparable to the net MHD effort in the calculation. In our experience, direct solution of eq. 
p-l| requires a momentum resolution pj + i/pj<^1.05, so we would need ~ 250 momentum bins 
(to say nothing of the spatial resolution issues emphasized earlier). Thus, the cost of following 
the electrons would exceed the cost of the MHD by at least an order of magnitude. Since these 
are computationally expensive simulations anyway, especially in 3D, that would make ordinary 
schemes completely impractical at present. 
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FIGURE CAPTIONS 



Fig. 1 Dynamical flow properties of the simulated MHD jet at t = 10.67. Normal gray scale 
is used, with high tones representing high values. Shown top to bottom are: (a) Log gas 
density, (b) Log magnetic pressure, (c) Gas compression, with shocks showing high tones. 

Fig. 2 (a,b): Log synchrotron emissivity computed at 1.4 GHz from the simulated electron 

distributions for Models 1 and 2 at the same time and with the same MHD properties as 
in Fig. 1. (c): "Pseudo emissivity", j c , using only MHD variables. The same dynamic range 
is shown for all; namely, 3 x 10 3 below the maximum value in each case. 

Fig. 3 Synchrotron emissivity spectral index distributions, ot\ and ot2, for the electron acceleration 
Models 1 and 2, as shown in Fig. 2. (a) Model 1, (b) Model 2. Regions are included only 
if they show in Fig. 2 and have a < 0.8. 

Fig. 4 (a): Log synchrotron emissivity at 1.4 GHz computed from Model 3, which includes 
strong synchrotron aging effects. The dynamic range is increased over Fig. 2 by a factor 
10 to 3 x 10 4 to capture the same physical regions, (b): Synchrotron emissivity spectral 
index distribution, 03 for Model 3. Regions are included only if they show in Fig. 4a. (c): 
Synchrotron spectral curvature between 1.4 GHz and 5 GHz for the same regions shown in 
Fig. 4b. 

Fig. 5 Scatter plot displaying the distribution of spectral index vs emissivity [(a), (c), (e)] and 
spectral magnetic field vs index [(b), (d), (f)] values for Models 1, 2 and 3 at t = 10.67 
Each point represents one zone on the simulation grid, using a stride of 2 to reduce the size 
of the plot file. The number of points in a given range on these plots can be seen as the 
area on the computational grid with those properties. The magnetic field is normalized by 

-^background • 

Fig. 6 A 2-D test calculation of the relativistic electron transport scheme described in section 
A. 3. It involves the evolution of two obliquely intersecting shocks, as described fully in the 
text. The computational box was given dimensions 1.5 x 1.5 and the time shown is t = 0.05. 

(a) Gas density, with the five dynamical flow regions labeled. High tones are high densities. 

(b) Electron momentum spectral index. High tones are flatter spectra. 

Fig. 7 Tests of the electron transport scheme based on simulations identical to Model 3 (see 
Table [l]), except these were done at half the spatial resolution. Gray scale image: Log of 
the relativistic electron density, (a),. . .,(d): Computed electron momentum distributions at 
indicated locations using 4 (triangle), 8 (star), and 12 (square) momentum bins to cover the 
range shown. The plot (b) was taken from a region subject only to synchrotron aging effects, 



and shows the results when t = t so (see eq. [3-11, with p = 10 4 , along with the analytic 
solution (solid curve). 
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